% This file take the steady-state the file "todynare_Reiter.mat"
% It writes the dynare file "code_dynare_rules.mod" which is the dynamic model
% to use perturbation techniques
% It doubles check that the steady state is OK, thanks to the "resid" function

% Be careful : We hard wire the fact fact there are 3 ex-ante different
% agents (only to aggregate their consumption).


clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

isOctave && warning('off', 'Octave:array-as-logical');


load ../todynare_Reiter.mat

aGrid = kron(ones(eco.ny,1),eco.aGrid);

% Constructing relevant objects
Ntot  = eco.na*eco.ny;    %total number of bins
ytype = kron(eco.ys,ones(eco.na,1)); % productivity type by bin 
yind  = int8(kron(double([1:eco.ny]),ones(eco.na,1))); % productivity index by bin 
na = eco.na;
ny = eco.ny;
D_ss = solution.stationaryDist(:);



K     = solution.K(1);
Ltot  = solution.L(1);
alpha = eco.alpha(1);
delta = eco.delta(1);
Y     = K^alpha*Ltot^(1-alpha);
FL    = (1-alpha)*K^alpha*Ltot^(-alpha);


coefB = 0.1;


ys    = eco.ys;


rho_u = 0.1;
G0 = .01*(1-rho_u); 

l = (eco.chi*eco.tau*solution.w)^(1/(1.0/eco.phi + 1.0 - eco.tau));
taut =  (( 1.0/eco.phi + 1.0)*eco.tau)/( 1.0/eco.phi +1.0 - eco.tau);
cv = solution.gc(:);
lv = solution.gl(:);
av = solution.ga(:);

x = cv - (1/(eco.chi*(1/eco.phi + 1)))*(lv).^(1+1/eco.phi);



N1 = na*ny/3;
N2 = 2*na*ny/3;

C1 = sum(D_ss(1:N1).*cv(1:N1));  % consumption of each type of agent
C2 = sum(D_ss(N1+1:N2).*cv(N1+1:N2));
C3 = sum(D_ss(N2+1:Ntot).*cv(N2+1:Ntot));

P     = ones(Ntot,1);
P(find(av<1e-8)) = 0;


file = 'code_dynare_reiter_rules.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     variables and params     %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%% Note that we order the variables to simplify the computation of the
    %%% Theil index of inequality.
    
str = ['var\n'] ; fprintf(fid, str);


for h = 1:Ntot
    str = ['S',num2str(h),' '] ; fprintf(fid, str); 
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;
end;


for h = 1:Ntot
    str = ['a',num2str(h),' '] ; fprintf(fid,str);     %  end-of-period
   
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;

end;

str = ['u r w FK FL A B K TFP Ltot Ctot C1 C2 C3 C1t C2t C3t Y I It\n'] ; fprintf(fid, str);
str = ['ut Kt Ct Bt wt rt taut tauk Zt Yt Lt tkt kappa BoYt G l T\n'] ; fprintf(fid, str);

%str = ['r w   mu   \n'] ; fprintf(fid, str);
%str = ['tau     \n'] ; fprintf(fid, str);


for h = 1:Ntot
   % str = ['a',num2str(h),' '] ; fprintf(fid,str);     %  end-of-period
    %str = ['at',num2str(h),' '] ; fprintf(fid, str);   % beginning-of-period
    str = ['x',num2str(h),' '] ; fprintf(fid, str);
    %str = ['S',num2str(h),' '] ; fprintf(fid, str);
    str = ['ww',num2str(h),' '] ; fprintf(fid, str);
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;

end;

str = ['; \n\n'] ; fprintf(fid, str);
str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi abar delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        initialisation        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(eco.beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['abar','    = ',num2str(eco.a_min,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(solution.G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(eco.chi,formatSpec),';\n']; fprintf(fid, str);

%%
equa = 1; % Numbering equations
tol = 0*10^-8;  % threshold for transition, to avoid to consider very small transitions

str = ['\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       model definition       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['model;\n\n']; fprintf(fid, str);
for h = 1:Ntot
    
    %%%%%%%% Budget constraint
    str = ['x',num2str(h),' = (1/(chi*taut))*l^(1+1/phi)*(',num2str(ytype(h),formatSpec),')^taut + (1+r)*',num2str(aGrid(h),formatSpec),' + T -a', num2str(h), ...
           '; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
end
    
    
%%%%%%%%%%%%%%%%  WEALTH TRANSITION  %%%%%%%%%%%%%%
for i = 1:Ntot
     str = ['a',num2str(i),' =  ww',num2str(i),'*',num2str(aGrid(Vind(i)),formatSpec),'+ (1 -ww',num2str(i),')*',num2str(aGrid(Vind(i)+1),formatSpec),';']; fprintf(fid, str);
     str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

end;
str = ['\n\n']; fprintf(fid, str);

  
for h = 1:Ntot    
 %%%%%%%% Euler equation

if P(h)==1
    yi = yind(h); % current productivity status
    threse = 0;
          str = ['(x',num2str(h),')^-gamma  = beta*(1+r(+1))*('];fprintf(fid, str);
          for j = 1:eco.ny % next period productivity status
                  if (eco.Piy(yi,j)>threse) % does agent i reach agent j ? effective transition toward the bin j
                     str = ['+',num2str(eco.Piy(yi,j),formatSpec),'*( ww',num2str(h),'*x',num2str(Vind(h)+(j-1)*na),'(+1)+(1-ww',num2str(h),')*x',num2str(Vind(h)+1+(j-1)*na),'(+1) )^-gamma']; fprintf(fid, str);
                  end;
          end;
            str = [')+',num2str(resE(h),formatSpec),';']; fprintf(fid, str);
            str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
else
        str = ['a',num2str(h),' =  ',num2str(av(h),formatSpec),';']; fprintf(fid, str);str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
end

end



%%%%%%%%%%%%%%%%%%%%%%%%%  evolutions of shares  %%%%%%%%%%%%%%%%
 str    = [' %%%%%%%% equation ',num2str(equa),': Shares \n\n']; fprintf(fid, str);
for i = 1:na
    for j=1:ny
        pass=0;
        ind =   i + (j-1)*na;
        str    = ['S',num2str(ind),' = ']; fprintf(fid, str);

        for ip = 1:na
            for jp=1:ny
             indp =   ip + (jp-1)*na ;
               if (Vind(indp)==i)&(eco.Piy(jp,j)>0) % goes from indp to ip (i+1 appears later)
                        str    = ['+S',num2str(indp),'(-1)*',num2str(eco.Piy(jp,j),formatSpec),'*(',num2str(aGrid(Vind(indp)+1),formatSpec) ,'-a',num2str(indp),'(-1))/',num2str(aGrid(Vind(indp)+1)-aGrid(Vind(indp)),formatSpec)]; fprintf(fid, str);pass=1;
               end
               if (Vind(indp)==i-1)&(eco.Piy(jp,j)>0) % goes from indp to ip (i+1 appears later)
                    if aGrid(Vind(indp))>0
                        str    = ['+S',num2str(indp),'(-1)*',num2str(eco.Piy(jp,j),formatSpec),'*(a',num2str(indp),'(-1)-',num2str(aGrid(Vind(indp)),formatSpec) ,')/',num2str(aGrid(Vind(indp)+1)-aGrid(Vind(indp)),formatSpec)]; fprintf(fid, str);pass=1;
                    else
                        str    = ['+S',num2str(indp),'(-1)*',num2str(eco.Piy(jp,j),formatSpec),'*(a',num2str(indp),'(-1)' ,')/',num2str(aGrid(Vind(indp)+1)-aGrid(Vind(indp)),formatSpec)]; fprintf(fid, str);pass=1;
                    end
              end
            end
       end
         if pass==0 % not assigned (constant share)
          str    = [num2str(D_ss(ind),formatSpec)]; fprintf(fid, str);
        end
       str    = [';']; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    end
end




%%%%%%%%%%%%%%%% Aggregate Constraints %%%%%%%%%%%%%%%%

%%%%%%%% Resource contraint
str = ['G = Gss*(1+u);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['I = K - (1-delta)*K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['G + Ctot + I = Y;\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%str = ['G + (1+r)*B(-1) + (Ctot +A - (1+r)*A(-1))  = Y +B - (r+delta)* K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Financial market clearing
str = ['A = B + K; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Production side
str = ['FL = (1 - alpha)*TFP*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*TFP*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['TFP = 1; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = TFP*K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate labor supply
str = ['Ltot = l*('] ; fprintf(fid, str);
for j = 1:ny
    str = ['+ ', num2str(eco.Sy(j),formatSpec),'*',num2str(eco.ys(j),formatSpec),'^(1+ taut/(1/phi+1))'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [');   %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;



%%%%%%%% Aggregate consumption   %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
str = ['Ctot = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+ S',num2str(h),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(ytype(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['C1 = 0'] ; fprintf(fid, str);
for h = 1:N1
    str = ['+ S',num2str(h),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(ytype(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['C2 = 0'] ; fprintf(fid, str);
for h = (N1+1):N2
    str = ['+ S',num2str(h),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(ytype(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['C3 = 0'] ; fprintf(fid, str);
for h = (N2+1):Ntot
    str = ['+ S',num2str(h),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(ytype(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Wage
str = ['w =  (1/chi)*(1/taut +1/(1/phi+1) ) *l^( (1+1/phi)*(1+1/phi)/(1/phi+1 +taut) )  ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Taxes
str = ['tauk = 1-r/FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['kappa = w/(FL^( (1/phi+1)*taut/(1/phi+1+taut) ));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate savings
str = ['A = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+ S', num2str(h),'*a',num2str(h)] ; fprintf(fid, str);
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%%% Fiscal system

str = ['tauk  = steady_state(tauk);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['taut = ',num2str(taut,formatSpec),';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

str = ['kappa = steady_state(kappa) -',num2str(coefB,formatSpec),'*(B - steady_state(B)) ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['T = 0 ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%%%%%%%%%% log-deviation variables %%%%%%%%%%%%%%%%
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(solution.K(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(solution.C(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['C1t = 100*(C1/',num2str(C1,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['C2t = 100*(C2/',num2str(C2,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['C3t = 100*(C3/',num2str(C3,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Bt = 100*(B/',num2str(solution.B(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(solution.w(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - ',num2str(solution.R(1)-1,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['tkt = 100*(tauk -',num2str(eco.tauk,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Zt = -ut;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(solution.L,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(solution.Y(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['It = 100*(I/',num2str(delta*K,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['BoYt = 100*(B/Y*',num2str(solution.Y/((1+eco.tauc)*solution.A-solution.K),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n'] ; fprintf(fid, str);
%%%%%%%%%%%%%%%% end of the model part


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



str = ['steady_state_model;\n\n'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['a',num2str(h),'        = ',num2str(av(h),formatSpec),';\n'] ; fprintf(fid, str);   
    str = ['x',num2str(h),'        = ',num2str(x(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['S',num2str(h), ' = ',num2str(D_ss(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['ww',num2str(h), ' = ',num2str(Wind(h),formatSpec),';\n'] ; fprintf(fid, str);
    
end;


str = ['l     = ', num2str(l,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = ', num2str(solution.R-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(solution.w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(1/eco.beta-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ltot  = ', num2str(solution.L(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['TFP   = ', num2str(1),';\n'] ; fprintf(fid, str);
str = ['A     = ', num2str(solution.A(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['B     = ', num2str(solution.A(1)-solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['tauk  = ', num2str(eco.tauk,formatSpec),';\n'] ; fprintf(fid, str);
str = ['taut   = ', num2str(taut,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(solution.C,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(solution.Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['kappa = ', num2str(eco.kappa,formatSpec),';\n'] ; fprintf(fid, str);
str = ['I  = delta*K;\n'] ; fprintf(fid, str);

str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['C1t    = 0;\n'] ; fprintf(fid, str);
str = ['C2t    = 0;\n'] ; fprintf(fid, str);
str = ['C3t    = 0;\n'] ; fprintf(fid, str);
str = ['Bt    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Tt    = 0;\n'] ; fprintf(fid, str);
str = ['tkt   = 0;\n'] ; fprintf(fid, str);
str = ['Zt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['It    = 0;\n'] ; fprintf(fid, str);
str = ['BoYt  = 0;\n'] ; fprintf(fid, str);
str = ['T  = 0;\n'] ; fprintf(fid, str);

str = ['G     = Gss;\n'] ; fprintf(fid, str);

str = ['C1  = ', num2str(C1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['C2  = ', num2str(C2,formatSpec),';\n'] ; fprintf(fid, str);
str = ['C3  = ', num2str(C3,formatSpec),';\n'] ; fprintf(fid, str);



str = ['end;\n\n'] ; fprintf(fid, str);
%%%%%%%%%%% end of the steady state part

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);
 

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       Aggregate shock        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 str = ['shocks;\n var eps =',num2str(G0,formatSpec),';\nend;\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%         Stoch simul          %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   

str = ['stoch_simul (order=1,irf=500,periods=1000) ut Bt taut kappa tauk Yt Ct Lt C1t C2t C3t It T;'] ; fprintf(fid, str);

isOctave && fflush(fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        Launch Dynare         %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    command = ['dynare ',file, ' noclearall'];
    eval(command);

name = ['fig_Reiter',num2str(rho_u),'.mat'];

save (name,'C1','C2','C3','ut_eps','Bt_eps','Yt_eps','It_eps','Lt_eps','Ct_eps','C1t_eps','C2t_eps','C3t_eps','kappa_eps','oo_')






